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We propose a description of nonequilibrium systems via a simple protocol that combines exchange- 
correlation potentials from density functional theory with self-energies of many-body perturbation 
theory. The approach, aimed to avoid double counting of interactions, is tested against exact results 
in Hubbard-type systems, with respect to interaction strength, perturbation speed and inhomo¬ 
geneity, and system dimensionality and size. In many regimes, we find significant improvement over 
adiabatic time dependent density functional theory or second Born nonequilibrium Greens function 
approximations. We briefly discuss the reasons for the residual discrepancies, and directions for 
future work. 
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Hybrid methods are a valuable option in physics, to 
merge concepts and perspectives into a more general 
and effective level of description. This work adds an 
item from condensed matter physics to the list; we pro¬ 
pose a hybrid method which combines non-perturbative 
exchange-correlation (XC) potentials from Time Depen¬ 
dent Density Functional Theory (TDDFT) [IHl] with 
many-body perturbative self-energy schemes from Non- 
Equilibrium Green’s Functions (NEGF) [IH^, to deal 
with systems with strong electronic correlations and out 
of equilibrium. 

An accurate first-principles description of the real¬ 
time dynamics of systems with strong electronic corre¬ 
lations is an important, difficult and basically unsolved 
problem of condensed matter research. General frame¬ 
works like TDDFT and NEGF do indeed allow for an 
in-principle-exact treatment of strong electronic correla¬ 
tions. However, they both rely on key ingredients [the 
exchange-correlation (XC) potential for TDDFT and the 
self-energy E for the NEGF] that in general are only 
approximately known. For TDDFT, a systematic and 
general way to include non-local, non-adiabatic effects 
in the XC potential is lacking, while for NEGF a main 
hindrance is that self-energies based on many-body per¬ 
turbation theory, already computationally demanding, 
are usually inadequate for strong electronic correlations. 
While considerable progress has been made for model sys¬ 
tem far away from equilibrium (see e.g. [51114)) or for the 
ah initio description of near-equilibrium situations (see 
e.g. mi El), a reliable first-principles treatment of the 
far-from-equilibrium regime is still lacking. 

Here, we suggest a step towards the solution of 
this problem, by a novel combination of TDDFT and 
NEGF, where perturbative (but systematic) memory- 
effect corrections augment a non-perturbative local adia¬ 
batic treatment of electronic correlations. The approach 
is fully conserving in the Kadanoff-Baym sense m and, 
using the so-called generalized Kadanoff-Baym ansatz 


[15] (see below), can be made viable for realistic systems. 

Putting in practice our proposal at the ab-initio level 
requires access to continuum non-perturbative XC poten¬ 
tials, and this point is addressed at the end of the paper. 
However, the scope of our method can already be illus¬ 
trated here using simple lattice models. This has the ad¬ 
vantage of avoiding complex implementations and tech¬ 
nicalities that, indispensable to deal with real-world sys¬ 
tems, are usually unnecessary (possibly even unwanted) 
for an explorative assessment of a new methodology. Our 
results show that in many situations (see also the sup¬ 
plementary material, SM) the hybrid method provides 
significant progress over adiabatic-TDDFT and pertur¬ 
bative schemes for NEGF, thus holding promise for an 
improved treatment of the nonequilibrium dynamics of 
realistic correlated systems. 

Models systems.- We consider small Hubbard-type ID 
and 3D clusters, isolated or coupled to two ID semi¬ 
infinite non-interacting leads. In the latter case, the clus¬ 
ter consists of 1 site (single impurity). These systems 
are exposed to time-dependent (TD) local perturbations 
and/or (where applicable) to electric biases in the leads. 
The Hamiltonian for the above setups is 

H = H, + Hi+H,i, (1) 

which has contributions from the cluster, the leads, and 
the cluster-leads couplings. In standard notation, 

H 'y ] J- y ] Ciit^fl; T y ^ Gn2^7li^,(2) 

{ij)^C,(T i^C zGC 

where (ij) labels nearest-neighbour sites in the cluster 
C, V' > 0 is the tunneling amplitude, ei{t) are time- 
dependent on-site energies in the cluster, and Ui are 
contact-interaction strengths. Further, fii = fii^ -\- 
For the lead Hamiltonian, Hi = Ha, where a = R{L) 
refers to the right (left) lead, and 

Ha = -V Y. + ( 3 ) 

{ij)^(y,(T O' 
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FIG. 1. XC potentials from the ID- (left) and 3D- (right) 
homogeneous Hubbard model. 

Here, ba{t) is the (site-independent) bias in lead a,V > Q 
the tunnelling amplitude and N^, = The cou¬ 

pling between the leads and the cluster (impurity) are 
given by 

Hcl = —V[ink. r -|- ^c ,<y) + (4) 

(7 

All energies units are expressed in terms of the hopping 
parameter V' (for the 1-site impurity cluster we use V 
instead), and time is measured in the units of the in¬ 
verse hopping parameter (assuming atomic units). We 
now switch to continuum variables for generality and 
notational convenience, and provide some elements of 
TDDFT and NEGF relevant to our approach. 

NEGF.- The nonequilibrium propagator G(l,2) = 
G(<i,ri,t 2 ,r 2 ) satisfies the equation of motion [i9t^ — 
/i(l)]G(l, 2) = (5(1, 2) -I- E(l, 3)G(3, 2)d3 (and a similar 
one for 12 )- Here, h = T + vh + Vext is the single-particle 
Hamiltonian, with kinetic energy T, Hartree potential 
vh, and external potential Vext- S = Eemti + Sa;c[G] is 
the self-energy, which introduces a memory dependence. 
We integrate over the Keldysh contour 7 II IS]- ^emb is 
an embedding self-energy which accounts for the leads 
(if present), while 'Exc accounts for XC effects [20]. Stan¬ 
dard approximations for Exc are Second-Born (2BA), T- 
matrix (TMA) and screened interaction (GW) [SlIT]. For 
real time, the lesser part of G (denoted G<) gives the 
density n{t,r) = —iG<(t, r, t, r) and the current. 

TDDFT. - The time-dependent density uks is obtained 
in terms of the Kohn-Sham (KS) orbitals (/(^(t, r). These 
obey the KS equation [T-|-r;*:s(t, r)](^K(t, r) = i9t(^,j(t, r), 
where vks = vh + Vext + Vxc, and Vxc accounts for XC 
effects. Then, nigs(t,r) = Within a 

NEGF treatment, the KS density can be obtained from 
[ii9ti-/i(l)-U2;c(l)]Gifs(l,2) = (5(1,2), with n^s(t,r) = 
—iG^g(t, r, t, r). In practical implementations, the func¬ 
tional dependence of Vxc on n is often replaced by an 
Adiabatic Local Density Approximation (ALDA), i.e. 
Va:ci[n],r,t) « vZHnir,t)). 

A hybrid TDDFT-NEGF approach. - Our proposal is to 
augment a perturbative self-energy E^^ from a conserv¬ 
ing many-body scheme with a non-perturbative XC po¬ 
tential local in space/time. Alternatively, this pre¬ 
scription can be be seen as recasting an ALDA-TDDFT 


based on in a NEGF approach, but augmenting it 
with a non-local, non-adiabatic perturbative self-energy 
E^^. To avoid double counting we subtract an ALDA 
potential obtained from the same approximation as 
was used for E^^. The basic equation of our approach is 

[ idt, - Ml) - <Mn(l)) + ] G(l, 2) 

= 5(1,2)+ /'srj(l,3)G(3,2)(i3. (5) 
•I7 

To actually proceed with Eq.®, at t = 0 we solve 
for G, i.e. we find v'^^[G], v^^G] and E^J[G\ self- 
consistently on the imaginary-time track; then we propa¬ 
gate G self-consistently on the Keldysh contour, thus ful¬ 
filling the conservations laws of Kadanoff and Baym. The 
hybrid scheme involves no additional computational costs 
compared to standard NEGF time propagation. Since 
the augmentation — (t) is of the form of a time- 

local potential, our scheme can similarly be implemented 
in a density matrix formalism. This means that a Gen¬ 
eralized Kadanoff-Baym Ansatz (GKBA) [THUHIIIS] can 
be employed to reduce computational costs allowing for 
first-principles calculations of realistic systems. 

The non-perturbative XG potentials- For lattice sys¬ 
tems, depends on the system’s dimensionality. In ID, 
we describe the non-perturbative, adiabatic local correla¬ 
tions in terms of v^P{t,r) « v^^^^^{n{t,r)) [21], and in 
3D in terms of v'^P(t,r) « {n{t,r)) [^. 

is computed with the Bethe-ansatz from the ID Hubbard 
model [MlEl], and with DMFT [28l|29] from the 

3D homogeneous Hubbard model ^5] . 

The correction.- For concreteness, in this paper 
E^^ and are computed in the 2BA (some results 
in the TMA are also shown). The calculation and use 
of for Hubbard-type interactions in a NEGF time 
evolution has been discussed before (see e.g. [50]) and is 
not repeated here. Rather, we provide additional details 
of the perturbative correction . For the homogeneous 
(Hubbard) reference system, we use (n) = , 

where (n) = (n) — To{n) — Enin), and the three 

terms on the RHS respectively are the total energy in 
the second Born approximation, the non-interacting ki¬ 
netic energy and the Hartree energy for the ID (or 3D) 
homogeneous Hubbard model. We compute Ef^ {n) in 
q) -space: 

Etl = f du} f dqImG^(w,q)/(w)(a; + eq), 

GTT) J-ao JbZ 

with G^ the retarded propagator, / the sta¬ 
tistical Fermi factor (we consider zero tempera¬ 
ture), Cq the single-particle energies, and n = 
( 2 ^ 0+1 I-ood^ Ibz dqlmG^{u}, q)/(w). In Fig.j^we plot 
for the ID and 3D Hubbard model, for different in¬ 
teraction values. We also show the non-perturbative po¬ 
tentials used in Eq. ([^. They exhibit 
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FIG. 2. Time dependent density at the central site of a 5®- 
site cluster (panel c) for U = 8 (top row) and t/ = 24 (middle 
row), for slow (left column) and fast regime (right column). 
The effective cluster is displayed in panel d. The interacting 
and perturbed site of the cluster is coloured in orange. The 
perturbations Vs^t are shown in the bottom row. 



FIG. 3. Time evolution of the density m at site 3 in a L = 8 
Hubbard ring with U — under the perturbation Vext(l,t). 
The parameters for the perturbation profile are to = 0, cr = 1 
for the step (s), to = 5.5,0- = 0.5 for the ramp (r) and to = 
2.5,0- = '/OA for the gaussian (g). 


a discontinuity at half-filling, which is always present in 
ID but only for large U values in 3D, reflecting the Mott- 
Hubbard metal-insulator transition |25j . The discontinu¬ 
ity is absent in the 2BA. Note that, at exactly half-filling, 
and are both zero. Finally, from the TMA 
is shown. The discontinuity is absent also in this case, 
and at low/high filling approaches 

Closed systems: the 3D case. - We start our analy¬ 
sis with a 3D cubic cluster with 5^ sites, open bound¬ 
ary conditions, and a single interacting (and perturbed) 
site at the cluster centre (Fig. [^). We compare time- 
dependent densities from the hybrid-approach, 2BA and 
ALDA, against exact results. The system is highly inho¬ 
mogeneous, and despite the local character of the inter¬ 
action and external perturbation, non-local effects are 
important: the exact Vxc (not shown) can have large 
nonzero components at all sites [25]. Using symmetry, 
we map the cluster to a 10-site one (Fig. il) , as in pS] . 
We consider both weak (U = 8 , panels a,b) and strong 
correlations {U = 24, panels c,d). The temporal shape of 
the external fields we use is Gaussian (bottom-row pan¬ 
els, red curves), with a slower or faster onset/offset (in 
the following, referred to as fast or slow perturbations). 
For additional time profiles we refer to the SM. 

For the weakly correlated, slowly perturbed case (panel 
a), all approximations follow the exact solution. For the 
fast perturbation (panel b), non-adiabatic effects emerge, 
and this leads to the failure of the ALDA; the remain¬ 
ing approximations perform well, with the hybrid method 
marginally better than 2BA. In contrast, for the slow per¬ 
turbation and stronger correlations (panel c), the agree¬ 
ment of the 2BA is poor, while the other treatments still 
follow the exact solution. For the most unfavourable and 
extreme regime of strong correlations and fast perturba¬ 
tions (panel d), ALDA and 2BA are largely out of phase, 
and only the hybrid approximation reproduces the main 
structures of the exact solution with the correct phase. 


Overall, the hybrid approximation exhibits a fairly good 
agreement in all regimes, and is superior to the others in 
the most extreme regime. 

Closed systems: the ID case.- We next consider when 
all sites are interacting and exposed to a space- and time- 
dependent perturbation. A 3D system for this situation 
which is also an exactly solvable benchmark is not easily 
accessible, due to the unfavourable scaling of the con¬ 
figuration space. We thus turn to a numerically more 
convenient ID test-case (this also makes possible to as¬ 
sess the hybrid approach at low-dimensionality), choos¬ 
ing a ID ring with 8 interacting sites (Fig. |^. To ex¬ 
plore the role of space inhomogeneity, we resort to a 
(rather artificial) perturbation sinusoidally modulated in 
space: Vextil.t) = sin{^l (j)k) F(^^), where \k = 2 ^= 
[k = 1,2, ...,4) and F is temporal profile. The phase (pk 
guarantees that the sine nodes are between sites and the 
amplitude at site I = 3 has always the same sign. For 
the time profile, F{t) = 0{t) (step, s), F{t) = l/(l-|-e“‘) 
(ramp, r) or F{t) = exp{—t^) (gaussian, g). Results are 
shown in Fig.(for a more systematic study see the SM). 

With highly inhomogeneous fields (Ai,A 2 ) no approx¬ 
imation reproduces the exact dynamics. Moreover for 
the hybrid method shows artificial density oscilla¬ 
tions. The latter, also present in the the TDDFT-ALDA 
based on are induced by the sharp discontinuity 

in Vx^^^^ and are not removed by the 2BA self-energy 
(thus, non-local, non-adiabatic effects beyond the 2BA 
should be also taken into account). For more homoge¬ 
neous fields (A 3 ) the different approximations compare 
more favourably to the exact dynamics with superior¬ 
ity of the hybrid method. Looking at 5 ^= 3 , the hybrid 
approximation is in phase with the exact curve but, for 
densities changing across half-filling, it still exhibits the 
artificial oscillations (see the SM). Further, ALDA does 
not perform well, and 2BA tends to be out-of-phase with 
the exact solution. Finally, for a slowly varying-in-space 
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FIG. 4. Single-impurity, one-orbital Anderson model with 
U = 2 (shown in panel e). a): Linear conductance G in the 
wide-band-limit for Twbl ~ 0.09 (strong correlations). The 
exact G is displayed, together with the Hartree-Fock (HF), 
2BA, ALDA and hybrid-method results. The density/spin- 
channel n/2 at the impurity is also shown (dashed line), n/2 
and G share the same vertical scale (in different units), b-e): 
Time dependent density n (b, c) and average current = 

(d, e) for the Anderson impurity with constant impurity gate 
voltage Vgate = eo = 0.25 and bias bL{t) = 0.56(f). The 
hopping parameter in the leads is 1/ = 1, the impurity-lead 
coupling is Vii„k = 0.5 (b,d) and Viink = 0.3 (c,e). 

perturbation (A 4 ) the hybrid approach (in contrast to 
the other approximations) is in excellent agreement with 
exact results. This applies for all time profiles g, s, r. 

Open systems - Finally, we test the hybrid method in 
open systems (Fig. |^. Specifically, using a single-orbital 
Anderson impurity coupled to two ID semi-infinite leads 
pTT] (system shown in Fig.|^), we consider i) the conduc¬ 
tance G in the wide-band limit (WBL), Fig. |^; ii) the 
finite-bias, finite-lead-width regime, Fig. &-e. Starting 
with i), we find the exact density (and thereby the exact 
linear conductance via the Friedel sum rule) in the WBL 
Fig. 1^ displays for [/ = 2 the absolute devi¬ 
ation A from the exact G as function of Vgate and for 
different approximate treatments. We consider stronger 
correlations {T^bl = = 0.09, see the plateau in 

the conductance); here except for 0.15 < n/2 < 0.28, 
the hybrid method performs as the best compared to 
2BA or ALDA, and it is significantly better in the range 
0.28 < n/2 < 0.42 (symmetrical considerations apply 
above half-filling). 

ii) Next, we consider ID tight-binding leads (of band¬ 
width 4F). We fix a static Vgate to to be away from 


the particle-hole symmetric ground state (where = 
= 0). As benchmark, we use open-ended, Anderson- 
impurity finite chains with up to L = 96 sites treated 
with tDMRG PU [37]. When = 0.5 (panels 

b,d), the agreement between hybrid and tDMRG den¬ 
sities/currents is fairly good, especially in the transients 
(n and j from tDMRG never fully reach a steady state 
within the simulation time, in contrast to hybrid, 2BA, 
and ALDA, ones However, for stronger correlations and 
lower transparency U/Vunk = 2/0.3 (panels c,e), the 
impurity density from the hybrid scheme is closest to 
the tDMRG one than other schemes, whilst for the cur¬ 
rents ALDA performs best. The unconvincing perfor¬ 
mance of the hybrid approximation for U/Vunk = 2/0.3 
comes probably from multiple-scattering processes, ne¬ 
glected by 2BA. To corroborate this conjecture we have 
tested the hybrid method also using the TMA, which in¬ 
cludes such processes. In Fig. and e) the TMA hybrid 
method shows an improvement over the ALDA and the 
pure TMA calculation and thus supports the conjecture 
(for an expanded discussion and additional results, see 
the SM). 

Conclusions and outlook.- By merging elements of 
TDDFT and NEGF, we proposed a simple, easy to im¬ 
plement, nonequilibrium scheme aimed to improve the 
treatment of local non-perturbative correlation effects 
and, at the same time, to incorporate non-local, non- 
adiabatic effects. Results from Hubbard-type systems 
are quite encouraging. Taking a mildly optimistic stand, 
we can argue that our approach extends the applicabil¬ 
ity of ALDA-TDDFT and NEGF based on perturbation 
theory, thus providing a way forward to merge (strong) 
correlations and memory effects in general. On the other 
hand, one can certainly envisage situations where non- 
perturbative and non-local correlations are very impor¬ 
tant, and this is where perhaps corrections beyond the 
2BA (e.g., GWA or TMA or mixed, or other) could be 
employed. We note that Hubbard-type systems usually 
are challenging benchmarks to perturbative approxima¬ 
tions such as 2BA, TMA or GWA. The latter generally 
perform much better for continuum systems with long- 
range interactions. Thus, we speculatively suggest that 
our hybrid method could perform even better for real¬ 
istic systems. This is where the real merits of our pro¬ 
posal could possibly be: Using continuum XC potentials 
tailored for strong correlations (obtained from e.g. the 
strictly correlated approach [35MD] . where the disconti¬ 
nuities in Vxc manifest in a different way) and simplifi¬ 
cations for perturbative self-energies (such as the GKBA 
[ 1111211123 ]), our approach would be a leeway to an im¬ 
proved first principle treatments of realistic systems in 
nonequilibrium when strong local electronic correlations 
and memory effects play a role. 

We wish to acknowledge M. Puig von Friesen for dis¬ 
cussions in the early stages of this work. 
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